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Variational Monte Carlo analysis of the Hubbard model with a confining potential: 
one-dimensional fermionic optical lattice systems 
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We investigate the one-dimensional Hubbard model with a confining potential, which may 
describe cold fermionic atoms trapped in an optical lattice. Combining the variational Monte 
Carlo simulations with the new stochastic reconfiguration scheme proposed by Sorella, we 
present an efficient method to systematically treat the ground state properties of the confined 
system with a site-dependent potential. By taking into account intersite correlations as well 
as site-dependent on-site correlations, we are able to describe the coexistence of the metallic 
and Mott insulating regions, which is consistent with other numerical results. Several possible 
improvements of the trial states are also addressed. 
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1. Introduction 

Strongly correlated particle systems have attracted 
much interest. Among them, recent advances in laser 
cooling techniques make it possible to confine cold atoms 
in an artificial lattice, the so-called optical lattice. 1 ' 2 By 
tuning the amplitude and the direction of the laser, the 
hopping integral between sites and the strength of the 
confining potential can be controlled. Furthermore, by 
means of the Feshbach resonance, 3 the interaction be- 
tween particles is also controlled, which thereby stimu- 
lates intensive experimental investigations on the effect 
of particle correlations in the optical lattice. 

Bosonic systems have been discussed intensively 4-8 
since the discovery of the supcrfluid-Mott insulator tran- 
sition in the confined system with rubidium ions. 1 The 
ground state properties of the system have theoreti- 
cally been studied in detail. 9-16 Recently a degenerate 
Fermi gas with potassium ions (also lithium) is also re- 
alized, 17-24 which stimulates further investigations on 
strongly correlated fermionic systems in the optical lat- 
tice. 25-27 

Theoretical analyses of the fermionic systems have 
been done by means of the numerical techniques such as 
the quantum Monte Carlo (QMC) method 28-30 and the 
density matrix renormalization group (DMRG). 31,32 It 
has been elucidated that the metallic regions coexist with 
the Mott insulating regions in the one-dimensional (ID) 
confined system with large interactions. 28, 29,31 However, 
these numerical methods face serious problems when 
they are applied to the 2D and 3D systems. The DMRG 
method enables us to obtain the precise results only for 
the ID systems, and the QMC method usually suffers 
from the minus sign problems for a large-cluster system. 
Therefore, it is highly desirable to find a powerful method 
to investigate the fermionic confined systems systemati- 
cally. 

One of the potential methods suitable for this pur- 
pose may be the variational Monte Carlo (VMC) calcu- 
lation, 33-36 which allows us to study not only ID but 
also 2D and 3D systems in the same framework. How- 



ever, we again encounter some problems when the VMC 
method is extended to the confined systems with a site- 
dependent potential. One of the most serious problems 
is that a large number of variational parameters should 
be treated in order to correctly describe the ground state 
properties. Therefore, it is necessary to introduce sophis- 
ticated techniques beyond the standard ones used so far 
in variational theory. 

In this paper, we present an efficient method based 
on the VMC simulations that is tractable and applicable 
to the optical lattice systems. To overcome the above- 
mentioned difficulty inherent in the confined systems, 
we make use of a stochastic reconfiguration with Hes- 
sian acceleration (SRH) scheme to minimize the ground 
state energy in a given parameter space. This method 
was recently proposed by Sorella 37 to discuss long-range 
electron correlations in the periodic Hubbard system. By 
taking the ID Hubbard model with harmonic confine- 
ment as a typical example, we will demonstrate that the 
new stochastic VMC simulations work quite well for the 
optical lattice system with a confining potential. 

This paper is organized as follows. In §2, we introduce 
the model Hamiltonian and briefly summarize the VMC 
method with the stochastic reconfiguration scheme. In 
§3, we apply the VMC simulations to discuss the effect of 
the particle correlations in the confined system, by using 
several distinct trial states. A summary and discussions 
are given in the final section. 

2. Model and Method 

We consider a strongly correlated fermionic system in 
the optical lattice. Here, we treat fermionic particles with 
s = 1/2 spin, for simplicity, although a potassium atom 
40 K, for example, has an s = 9/2 spin. 26 In this pa- 
per, we restrict our discussions to the ID model in order 
to clearly demonstrate the potentiality of our numerical 
method applicable to the confined systems with a site- 
dependent potential. The simplified model Hamiltonian 
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we consider here reads 28, 29 



H 



<i,j>a 



C i<j C 3° 



+ ^2 V l n ia 



U^n^mi, (!) 



where c ia (c^) is a creation (annihilation) operator at 



the ith site with spin cr(=|,J,), and 



,c.«x- t is 



the nearest neighbor hopping matrix, and U(> 0) the 
on-site repulsive interaction. Interacting fermions are as- 
sumed to be trapped by the harmonic potential Vi , where 
Vl/2-i = Vl/2 = at two sites around the center, 
Vo = Vl-x = V at the edges (L is the length of the 
system). In the following, we fix t = 1 as an energy unit. 

In the confined system, the local particle count de- 
pends on the site because the confining potential Vi is 
site-dependent. Namely, the particle density around the 
center is larger than that around the edge. This sug- 
gests that the introduction of the repulsive interaction 
changes the distribution of various local quantities such 
as the particle density and the double occupancy. There- 
fore, the careful treatment of the particle correlations is 
desirable to clarify the ground state properties. 

To this end, we make use of the VMC method. 33-36 
In comparison with other numerical methods, it has an 
advantage in extending the calculation to the 2D and 3D 
fermionic systems, as mentioned above. This method has 
successfully been applied to various correlated electron 
systems with periodic potentials. 38-48 

We introduce here a trial state with Jastrow-type pro- 
jection operators as 



10 



|*) =J\$), 



(2) 



where |$) is the Slater determinant for the non- 
interacting state and the Jastrow operator J incorpo- 
rates two-particle correlations, whose explicit form will 
be given in the next section. In order to include the ef- 
fects of the inhomogeneous potential, we should deal with 
a large number of site-dependent variational parameters. 
This is contrasted to the standard problems in condensed 
matter physics, 33- where only one or two parameters 
are sufficient to describe the ground state properties. Al- 
though those parameters are invaluable to obtain a more 
accurate solution in the confined systems, it becomes 
much more difficult to determine their optimized values. 

For optimizing the parameters, there exist some iter- 
ative methods. In ordinary iterative methods, the set of 
parameters x can be iteratively changed by x — > x + 7. 
Here 7 is determined by minimizing the energy 



AE = E(x + 7) - E(x) 
= l~f T B~f - f T 7, 



(3) 



where f and B denote the first energy derivative and the 
Hessian matrix. As a result, 7 is denoted as 



7 = 5 -1 f. 



(4) 



In order to optimize the parameters efficiently, it is 
important to evaluate f and B. The SRH scheme, which 
is based on the iterative methods, allows us to optimize 
enormous parameters even more efficiently 37 In compar- 
ison with other iterative techniques such as the Powell's 



method and the quasi-Newton method, it has an advan- 
tage in reaching the ground state in a few iterations. The 
important idea in this scheme is that f and B are evalu- 
ated by using the statistical fluctuations of each projec- 
tion operator, and in the framework of VMC, the efficient 
and rapid convergence is achieved. Furthermore, we can 
introduce an additional free parameter (3 in the matrix 
-B, which can accelerate the convergence. This parameter 
originates from the expansion coefficient of the trial wave 
function, and contributes to the efficiency and the accu- 
racy of the calculation. This scheme was originally intro- 
duced to analyze a trial wave function with longer-range 
density-density or spin-spin correlations in the periodic 
systems. 37 

In this paper, we extend the SRH scheme so as to treat 
multivariable trial wave functions including the effects of 
the inhomogeneous potential. When the site-dependent 
parameters are introduced, the Hessian matrix B should 
become positive semidefinite with some zero eigenvalues. 
This reflects the fact that no particle exists around edges 
in the confining potential. In this case, some parameters 
do not contribute to the energy and thus cannot be de- 
termined straightforwardly. Therefore, it may be difficult 
to obtain B~ x in eq. (4). To overcome this difficulty, we 
propose to make use of 

7 = B+f (5) 

instead of eq. (4) . Here B + is the Moore-Penrose pseudo- 
inverse matrix of B. This matrix is always defined 
and has similar properties to the inverse matrix. We 
will demonstrate below that this new stochastic scheme 
works well and thereby enables us to discuss how parti- 
cle correlations affect the ground state properties in the 
confined system. 

3. Results 

In this section, we perform the VMC simulations with 
several trial states to clarify the role of the on-site and 
intersite particle correlations, by which we confirm the 
potentiality of our method. Here, we deal with the system 



(L = 100) with iV T 



30 and V = 10, where N a is 



the number of particles with spin a. 

3. 1 Gutzwiller paramagnetic state 

First, we consider the Gutzwiller paramagnetic trial 
state, which incorporates the on-site particle correla- 
tions. The Gutzwiller trial state 50-52 is explicitly given 
as 



V G 



cxp 



(6) 
(7) 



where |<&o) is the Slater determinant for the nonintcr- 
acting paramagnetic ground state. Note that onii = 
0, . . . , L — 1) is a site-dependent variational parameter 
with a,i = (XL—i—i, which may describe the effect of local 
correlations in our inhomogeneous system. 

By optimizing L/2(= 50) variational parameters {at} 
by means of the SRH scheme, we calculate local quanti- 
ties for the confined system. Let us first look at how ef- 
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(a) Particle Density 



10 15 
iteration 



Fig. 1. Calculated energies in each iteration. Filled circles, 
squares and diamonds are the results obtained from three dis- 
tinct trial states (see text) when U = 6 and V = 10. 



ficicntly our scheme works for optimizing the variational 
parameters. In Fig. 1, we show the calculated energy as 
a function of the iteration when it starts from a certain 
initial state (\^q) and \*&s) arc defined in the following 
subsections). For each iteration, we perform the VMC 
calculation with large samplings (~ 500, 000). An impor- 
tant point is that the ground state almost converges in 
a few iterations, although the energy for each trial state 
slightly oscillates even after sufficient iterations. We find 
that this does not depend on the number of the Monte 
Carlo samplings for each iteration, which implies that 
statistical errors in our calculations mainly come from 
this oscillation. Anyway we find that the stochastic VMC 
scheme with eq. (5) works quite well for our confined sys- 
tem. 

Several physical quantities thus calculated are sum- 
marized in Fig. 2. It is found that the introduction of 
the interaction slightly changes the local particle density 
(rii), and decreases the local double occupation (n^n^). 
At the same time, spin fluctuations — (Sf)) 2 ) are 
enhanced and density fluctuations ((rij — (fii)) 2 } are sup- 
pressed. This implies that the repulsive interaction U 
excludes doubly occupied states at each lattice site. It 
should be noticed that nearest-neighbor spin correlations 
{(SfS*_ 1 + SfSf_ 1 )/2) are strongly enhanced in the two 
regions around the 30th and 70th sites. In those regions, 
the average of the particle density is almost unity and 
density (spin) fluctuations are suppressed (enhanced) 
strongly. Namely, the repulsive interaction U has a ten- 
dency to localize particles and enhance local spin fluctu- 
ations, as should be expected. 

In this way, the trial state with site-dependent 
Gutzwillcr factors can take into account the local parti- 
cle correlations to some extent. In fact, the ground state 
energy obtained from the Gutzwiller trial state eq. (6) is 
lower than that obtained by the ordinary Hartree-Fock 
approximation, as shown in Fig. 3. However, we could not 
reproduce the Mott insulating region that should be spa- 
tially extended in a finite region, as suggested by other 
groups. 28 ^ 29 ' 31 



(c) Spin 
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(d) Spin Fluctuation 
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(e) Density Fluctuatioi 
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Fig. 2. (a) particle density, (b) double occupation, (c) short-range 
spin correlation, (d) spin fluctuations, and (e) density fluctua- 
tions, which are obtained by the Gutzwiller wave function | SE'c?) . 



3.2 Intersite doublon-holon correlations 

To take into account the particle correlations more 
precisely, we introduce the intersite doublon-holon (DH) 
correlations. 35 ' 36 The corresponding trial state is given 
as, 



l*0> = VqPg\$q) 



Vq = exp 



(8) 
(9) 



where 



Qi = Di JJ(1 - H l+T ) + Hi IJ(1 - A+t) (10) 

T T 

with Di = n^riii, Hi = (1 — ?i i |)(l — ?i^) and r runs over 
all the nearest neighbors. The variational parameter a 1 
is assumed to be independent on the sites, for simplicity 
(see the discussion in 3.3). Note that the doublon-holon 
projection operator Vq takes into account the nearest- 
neighbor correlations between a doubly occupied site and 
an empty site, which may be important to stabilize the 
Mott insulating state. 

The ground state energy thus computed is shown in 
Fig. 3. It is found that the energies obtained from the 
trial states eqs. (6) and (8) are almost the same in the 
weak coupling region. On the other hand, in the strong 
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(a) Particle Density 



Fig. 3. The ground state energy as a function of the repulsive in- 
teraction U. Open squares, open circles, open triangles, filled 
squares, and filled circles are the results obtained from the 
Gutzwillcr(GW), paramagnetic DH+GW, magnetically ordered 
DH+GW, SC+GW and SC+DH+GW trial state (see text for 
each notation). The energy obtained by the Hartree-Fock (HF) 
approximation is shown by the dotted line, for comparison. Sta- 
tistical errors originating from stochastic oscillations are less 
than the size of symbol for each point. 



coupling region U > 4, we find that the ground state 
energy for the trial state (8) is lower than the other, im- 
plying that the DH correlations play an important role in 
the region. In fact, characteristic behavior appears in var- 
ious physical quantities beyond U — 4, as shown in Fig. 
4. It is seen that when the interaction is increased, the 
plateau regions appear in the curve of the particle density 
at (rii) = 1, where the double occupancy is suppressed 
strongly. Furthermore, we find that spin (density) fluctu- 
ations are considerably enhanced (suppressed). This im- 
plies that the commensurability (half-filling condition) is 
induced sclf-consistcntly due to particle correlations and 
thus the Mott insulating states are stabilized in finite 
regions in contrast to the results of the Gutzwiller trial 
state (see Fig. 2). 

The emergence of the Mott insulating state is also seen 
clearly in the two-point spin correlation (— '(SfSj) 
(i = 40 fixed) shown in Fig. 5. At U = 4, the local par- 
ticle density, (m) ~ 1.2, at the 40th site is larger than 
half filling, which means that the metallic state is real- 



(c) Spin 




(e)Density Fluctuation 



site 



Fig. 4. (a) particle density, (b) double occupation, (c) short range 
spin correlation, (d) spin fluctuations, and (e) density fluctua- 
tions, which are obtained from the DH+GW trial state, IvPq). 
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Fig. 5. Staggered long-range spin correlations ( — l) l+J SfS| be- 
tween the 40th site and the j'th site. Squares, circles and dia- 
monds are the results for the system with U = 4.0, 5.0 and 6.5. 



ized around there, and therefore the spin correlations are 
extended up to a few neighbor sites. On the other hand, 
in the large U case, we find that the spin correlation 
length becomes much longer, signaling the formation of 
the Mott insulating state when U > 6. 

Summarizing, by including on-site and intersite parti- 
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cle correlations, our variational theory properly repro- 
duces the coexistence phase, which is consistent with 
those of QMC and DMRG simulations. 28 ' 29 ' 31 Further- 
more, by comparing the results obtained from the two 
trial states eqs. (6) and (8), we can say that intersite DH 
correlations are essential to describe the coexistence of 
metallic and Mott insulating states in the confined sys- 
tem. 

3.3 Further improvements of the paramagnetic state 

Here, we wish to briefly comment on further possible 
improvements of the trial wave function for the paramag- 
netic state. We examine a trial state which incorporates 
the effect of intersite spin correlations. The correspond- 
ing trial state we consider is 



l*s) 
\*sq) 

Vs 



Vs\*g), 
Ps\*q)- 



exp 



(11) 

(12) 
(13) 



where the spin correlation operator Si is defined as, 



£ = \ (s*su 



(14) 



The ground state energy calculated for each state is 
shown in Fig. 3. As seen in the figure, we cannot find 
a drastic change in the energy (as well as the physi- 
cal quantities) although a lot of variational parameters 
for spin correlations are introduced. The results indicate 
that the nearest-neighbor S z -S z spin correlations are not 
so important in comparison with the DH correlations. We 
think that Si ■ S^+i spin correlations or longer-range spin 
correlations are necessary to improve the trial state. 

One may also take into account site- dependent DH cor- 
relations to improve the trial state as, 



\^Q') = exp 



(15) 



In general, it is difficult to find the ground-state pa- 
rameters for this wave function in the strong coupling 
regime even by means of the SRH scheme. Nevertheless, 
by employing smaller lattice systems, we have checked 
that the site-dependent DH correlations give similar re- 
sults to that of the uniform DH correlations. Further im- 
provements of the iterative technique should allow us to 
discuss the effects of the site-dependent DH correlations 
more precisely, which is now under consideration. 

3.4 Examination of the magnetically ordered state 

We discuss here the ground state properties by ex- 
ploiting a trial state including the magnetic order. It is 
known that the Hubbard model in the periodic lattice 
without confining potentials has a tendency to stabilize 
the long-range magnetic order when the particle density 
satisfies the commcnsurability, i.e. the condition of half 
filling. In the ID Hubbard model, however, the real long- 
range order is not realized 53-56 due to large quantum 
fluctuations, in contrast to the 2D square-lattice system 



where the magnetic order is stabilized. 57 Nevertheless, it 
is expected that the effect of enhanced spin correlations 
can be incorporated to some extent even in terms of a 
Hartree-Fock type approximation. 

In our inhomogeneous system, the spatially-modulated 
magnetic properties may show up. Therefore, to examine 
a trial state with magnetic order, we start with the unre- 
stricted Hartree-Fock (UHF) state, where we determine 
the spatially modulated particle density and magnetiza- 
tion at each site self-consistently. 58 

The trial state for VMC simulations is then written as 



VqV g \$af(M)), 



(16) 



where \§af(M)) is the non-interacting state obtained by 
the UHF approximation with an additional variational 
parameter M, which corresponds to the total magne- 
tization. Note that the SRH algorithm can be applied 
only to the Jastrow-type variational parameters. There- 
fore, by performing the SRH scheme with a fixed M, we 
optimize the variational parameters ({aj}, a'). We show 
the particle density (n, ) and the average of the staggered 
magnetization \x in Figs. 6 and 7. The latter quantity is 
defined as. 



/' 



(17) 



We find that the magnetization is almost zero and little 
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Fig. 6. Filled (open) circles represent the density distributions of 
particles with f (J.) spin and a solid line the local particle density, 
when (a) U = 1.0, (b) U = 2.0, (c) U = 3.0, and (d) U = 4.0. 



difference appears in the spin-dependent density distri- 
butions when U < 2. This suggests that the paramag- 
netic ground state is realized in the case. On the other 
hand, the increase of the interaction yields the spatially- 
modulated magnetization, as shown in Fig. 6. An im- 
portant point is that the large staggered magnetization 
is induced in the vicinity of the i = 30th and i = 70th 
sites, where the plateau appears in the curve of the parti- 
cle density, as shown in Fig. 6. This implies that magnetic 
correlations are strongly enhanced in the finite insulat- 
ing region, which is consistent with those discussed in 
the previous subsection. 
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Fig. 7. The average of the staggered magnetization fi as a func- 
tion of U. Filled and open circles represent the magnetization 
obtained by the UHF approximation and the VMC method with 
the trial state |*q F ). 



Note that the VMC simulations with the HF approx- 
imation cannot describe a paramagnetic quantum liquid 
state correctly in ID. Nevertheless, we can see that the 
magnetically ordered state is more stable than the para- 
magnetic state in the large U region within our VMC 
analysis (see also Fig. 3). Therefore, if we can properly 
incorporate quantum fluctuations that suppress the mag- 
netic order in the Hartrce-Fock state, we end up with 
an excellent trial state. Alternatively, starting from the 
paramagnetic state, we should incorporate spatially ex- 
tended spin correlations more precisely. These problems 
are now under consideration. 

4. Summary and Discussions 

We have investigated correlation effects in the f D Hub- 
bard model with harmonic confinement, which may be 
relevant to correlated fermionic systems trapped in the 
optical lattices. We have aimed at developing an efficient 
method based on the VMC simulations that is tractable 
and applicable to the optical lattice systems. One of the 
central problems to be resolved is how to treat a large 
number of site-dependent variational parameters. 

As a first step to incorporate the effect of particle cor- 
relations in such inhomogeneous systems, we have exam- 
ined a trial state with site-dependent Gutzwillcr factors 
(on-site correlations). By exploiting the stochastic re- 
configuration with Hessian acceleration scheme, we were 
able to treat more than 50 Gutzwiller variational param- 
eters, whereas the standard Gutzwiller approach deals 
with only one or two parameters. We have thus clar- 
ified how the introduction of the repulsive interaction 
changes the distribution of local quantities. It turns out 
that although site-dependent correlations are taken into 
account, the Mott insulating region cannot be properly 
described only in terms of the on-site Gutzwillcr factors. 

We have thus improved the trial state by taking into 
account the intersite DH correlations. It has been eluci- 
dated that the coexistence state with metallic and Mott 
insulating regions emerges for large U in this case, which 
is consistent with those obtained by other numerical 
studies. This suggests that the intersite DH correlations 



are essential to describe the metal-insulator coexistence 
phase. We have further examined how the intersite spin 
correlations modify the results, and have found that as 
far as Ising-type spin correlations are concerned, the re- 
sults are little changed. Nevertheless, there still remains 
a possibility to improve the results by incorporating the 
Heisenberg-type isotropic spin correlations. It has been 
also shown that the variational treatment starting from 
the mean-field antifcrromagnctic state reasonably repro- 
duces the coexistence of metallic and insulating regions, 
although the real long-range order should be prohibited 
due to large quantum fluctuations inherent in ID sys- 
tems. In this way, the present scheme is quite flexible, 
thereby providing us with various ways to improve our 
variational treatment of the Hubbard model with a site- 
dependent potential. 

In conclusion, we have confirmed that the VMC sim- 
ulations with numbers of site-dependent parameters can 
be efficiently performed with the aid of the stochastic ac- 
celeration scheme. This demonstrates that the stochas- 
tic VMC method provides a powerful tool to treat 
the correlation effects in optical lattice systems with 
a site-dependent potential, which have not been easily 
tractable by conventional optimization procedures. This 
naturally motivates us to extend our discussions to the 
2D and 3D optical lattice systems, which is now under 
consideration. 
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